### Load library

library('SDMTools')

### Establish the out directory

out.dir = '/home1/99/jc152199/ChapterOne/AggregateFreqDistributions/'

### List the two spatial layers

accu = read.asc.gz('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/microCLIM_05.asc.gz')
anu = read.asc.gz('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/BC6_monthly/bc05.asc.gz')

### Bring in species occurrences

occurs = read.csv('/home1/99/jc152199/MAXENT/occurs/BIOCLIM_Occurrence_SWD.csv',header=T)

occurs = occurs[,c(1:3)]

### Read in raw physiology data from base.dir

minmax = read.csv('/home1/99/jc152199/MAXENT/output/PhysData/minmax.csv', header=T)

### Plot

for (spp in c('COPAENI','COPMONT','COPHOSM','COPINFA','COPBOMB','COPNEGL'))

	{
	
	### First obtain the phys results for the species of interest
	
	spp.mm = minmax[which(minmax$Species==spp),]
	
	### Now subset occurs to species of interest
	
	t.occur = occurs[which(occurs[,1]==spp),]
	
	### Extract climate from these occurrence points
	
	t.accu = extract.data(cbind(t.occur[,3],t.occur[,2]),accu)
	t.anu = extract.data(cbind(t.occur[,3],t.occur[,2]),anu)
	
	### Calculate density objects
	
	accu.density = density(t.accu)
	anu.density = density(t.anu)
	
	### Calculate the x and y lims
	
	xlims = range(c(accu.density$x,anu.density$x,spp.mm[1,5]+.5))
	ylims = c(0,max(c(accu.density$y,anu.density$y)))
	
	### Calculate a polygon that represents CTMax 95% CI (based on normally distributed data)
			
	max.poly = data.frame(x=c(spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96),y=c(ylims[2],0,0,ylims[2]))
			
	#### Now open the .png device driver and create a file for each species (but all in the same folder)

	png(paste(out.dir,spp,'.png',sep=''),units='cm', height=18, width=18, res=600)
	
	### Configure the plot space, and plot the BRT density object
	
	plot(accu.density, col='blue', xlim=c(xlims[1]-6,xlims[2]), ylim=ylims, sub='', main=paste(spp),xlab='Max Temp of the W@armest Period', ylab='Density')
	
	### Add the AWAP density object
	
	points(anu.density, type='l', col='red')
	
	### Add line CTMax
	
	points(c(spp.mm[1,5],spp.mm[1,5]),c(0,ylims[2]),type='l',lty=3, col='#FF0000')
	
	### Add a legend
	
	legend('topleft',legend = c('accuCLIM Distribution','ANUCLIM Distribution','CTMax 95% CI'),text.col=c('blue','red','black'),lty=c(1,1,1),col=c('blue','red','#FF000025'),lwd=c(1,1,8), bty='n')
	
	### Plot the polygon if it exists	
		
	polygon(x=max.poly[,1],y=max.poly[,2], col='#FF000025', border=NA)	
	
	### Close the device driver
	
	dev.off()
	
	}
	
# Close loop
	
	
	
